Exercise 20

stalkie <- read.csv(url("http://www.zoology.ubc.ca/~schluter/WhitlockSchluter/wp-content/data/chapter11/chap11e2Stalkies.csv"))


m0 = 9
M <- mean(stalkie$eyespan)
n <- length(stalkie$eyespan)
sigma <- sd(stalkie$eyespan)
SE <- sigma / sqrt(n)
test.statistic <- (M - m0) / SE

p.value = 2 * pt(-abs(test.statistic),df=8)
print(p.value)

## [1] 0.1327105

x=seq(-4,4,by=0.01)
y = dt(x,df=8)
plot(x,y,type="l")

abline(v = test.statistic,col="dodgerblue")
text(-3,0.2,"test statistic",col="dodgerblue")
text(2,0.3,"t-distribution, 
     with 8 degrees of freedom")

x1 = seq(-4,test.statistic,by=0.01)
y1 = dt(x1,df=8)

x2 = rev(x1)
y2 = y1*0

polygon(c(x1,x2),c(y1,y2),col = "red")
polygon(-c(x1,x2),c(y1,y2),col = "red")
text(-3.5,0.05,"P value",col="red")

The p value of the t-test is

test = t.test(stalkie$eyespan-9)

print(test$p.value)

## [1] 0.1327105

and the test statistic is

print(test$statistic)

##         t 
## -1.673772

We can summarize this in a table:

dat = data.frame(pValues = c(p.value,test$p.value),tStatistics = c(test.statistic,test$statistic))
colnames(dat) = c("P Value","Test Statistic")
rownames(dat) = c("by hand","using t.test()")
kableExtra::kable(dat) 

Exercise 21

data = c(58.9, 7.8, 108.6, 44.8, 11.1, 19.2, 61.9, 30.5, 12.7, 35.8, 7.4, 39.3, 24.0, 62.1, 24.3, 55.3, 32.7, 65.3, -19.3, 7.6, -5.2, -2.1, 31.0, 69.0, 88.6, 39.5, 20.7, 89.0, 69.0, 64.9, 64.8)

1. What is the sample size n?

n = length(data)

print(paste("The sample size is ",n))

## [1] "The sample size is  31"

2. What is the mean of these data points? Remember to give the units.

m = mean(data)

print(paste("The mean shift is ",m, "meters."))

## [1] "The mean shift is  39.3290322580645 meters."

3. What is the standard deviation of elevational range shift? (Give the units as well.)

sd = sd(data)

print(paste("The staandard deviation of the shift is ",round(sd,digits=2)," meters."))

## [1] "The staandard deviation of the shift is  30.66  meters."

4. What is the standard error of the mean for these data?

SE = sd/sqrt(n)

print(paste("The standard error of the mean is ",round(SE,digits=2),"."))

## [1] "The standard error of the mean is  5.51 ."

5. How many degrees of freedom will be associated with a confidence interval and a one-sample t-test for the mean elevation shift?

There are 30 (= \(n-1\)) degrees of freedom.

We can check with the t.test() function:

tt = t.test(data)
tt

## 
##  One Sample t-test
## 
## data:  data
## t = 7.1413, df = 30, p-value = 6.057e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  28.08171 50.57635
## sample estimates:
## mean of x 
##  39.32903

6. Calculate the 95 % confidence interval for the mean using these data.

We can either use our results from the t-test before:

tt$conf.int

## [1] 28.08171 50.57635
## attr(,"conf.level")
## [1] 0.95

or

upper = m + 2*SE
lower = m - 2*SE

c(lower,upper)

## [1] 28.31452 50.34355

A more precise version:

upper = m + qt(0.975,df = n - 1)*SE
lower = m + qt(0.025,df = n - 1)*SE
c(lower,upper)

## [1] 28.08171 50.57635

7. For the one-sample t-test, write down the appropriate null and alternative hypotheses.

\(H_0:\) The mean elevational shift is 0 meters. \(H_A:\) The mean elevational shift is not 0 meters.

8. Calculate the test statistic t for this test.

We canu again just use t.test():

tt$statistic

##        t 
## 7.141309

or do it by hand:

mu.0 = 0 # as specified in the null hypothesis
t.stat = (m - mu.0)/SE
t.stat 

## [1] 7.141309

1. What assumptions are necessary to do a one-sample t-test?

The data should be approximately normally distributed. We can check this with a histogram

hist(data)

A slightly more precise way would be to do a QQ plot

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

qqPlot(data)

## [1]  3 19

As long as the points are within the dashed lines, you can safely do a t-test.

Exercise 22

10. Calculate the P-value for this test as accurately as you can.

The most accurate way is either to use t-test again

tt$p.value

## [1] 6.056689e-08

We can also do it by hand. For this it is good to make a sketch first:

# we first calcualte the null distirbution of our t statisitc
x = seq(-8,8,by=0.01)
y = dt(x,df = n-1)

# and plot it
plot(x,y,type="l",xlab="t",ylab="density")

# we add the observed t statistic
abline(v = t.stat,col="slategray")

We see in the figure that the observed t statistic is far away from the mode of the t distribution. The p value shoul dtherfroe be very small. To calcualte it, we calcualte the area under the curve to the right of the statistic, and double this value to account for extreme events on the left end of the distirbution (we perform a two-sided test). The function pt would give us the area to the left of the test statistic, so we use 2*(1-pt(t.stat,df = n-1)) to get the p-value:

2*(1-pt(t.stat,df = n-1))

## [1] 6.056689e-08

11. Did species change their highest elevation on average?

We choose a signifcance threshold of 5 percent, i.e., we want to be 95 percent certain that we do not have a false positive. Th p-value is clearly smaller than that. Thus, the t-test shows that there is a statistically significant change in the highest elevation. Because the p-value is so small, we can be very certain of this result as it would alos hold at much lower signifcance thresholds.

Exercise 23

Previously we used a one sample t-test to calculate the confidence intervals. This test also tested the null hypothesis that the means of each group are 0. If we want to compare means between groups we need to do a two sample t-test:

t.test(stickleback$plates[stickleback$genotype=="mm"],
       stickleback$plates[stickleback$genotype=="Mm"])

## 
##  Welch Two Sample t-test
## 
## data:  stickleback$plates[stickleback$genotype == "mm"] and stickleback$plates[stickleback$genotype == "Mm"]
## t = -32.001, df = 208.06, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -41.09355 -36.32416
## sample estimates:
## mean of x mean of y 
##  11.67045  50.37931

t.test(stickleback$plates[stickleback$genotype=="mm"],
       stickleback$plates[stickleback$genotype=="MM"])

## 
##  Welch Two Sample t-test
## 
## data:  stickleback$plates[stickleback$genotype == "mm"] and stickleback$plates[stickleback$genotype == "MM"]
## t = -95.49, df = 167.89, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -52.16670 -50.05336
## sample estimates:
## mean of x mean of y 
##  11.67045  62.78049

t.test(stickleback$plates[stickleback$genotype=="Mm"],
       stickleback$plates[stickleback$genotype=="MM"])

## 
##  Welch Two Sample t-test
## 
## data:  stickleback$plates[stickleback$genotype == "Mm"] and stickleback$plates[stickleback$genotype == "MM"]
## t = -10.262, df = 207.1, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -14.78364 -10.01871
## sample estimates:
## mean of x mean of y 
##  50.37931  62.78049

We find that all t-test yield a significant result (small p values). This means that the null hypothesis of the mean of two groups being the same can be rejected for all comparisons.

A shorter way is by using the function pairwise.t.tests. The First argument is the variable on which the t-test should be performed. The second parameters is the grouping factor: the categorical variable that separates the first argument into different groups.

 pairwise.t.test(stickleback$plates,stickleback$genotype)

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  stickleback$plates and stickleback$genotype 
## 
##    mm      Mm     
## Mm < 2e-16 -      
## MM < 2e-16 1.5e-15
## 
## P value adjustment method: holm

How does this finding relate to the confidence intervals? The confidence intervals are not overlapping, strongly suggesting that the means between groups are indeed different (at the 95 percent level). Thus, the t-test and the confidence interval indicate the same result.

Exercise 24

mouse_cooling <- read.csv("~/Dropbox/Teaching/StatisticsForBiology/mouse_cooling.csv")

mouse_cooling = mutate(mouse_cooling,diff = freshly.killed - reheated)   
ggplot(mouse_cooling) + 
  geom_histogram(aes(x = diff),bins=7)

The data is clearly right skewed. This can also be seen with a QQ plot:

library(car)
qqPlot(mouse_cooling$diff)

## [1] 15  2

The null hypothesis is

\(H_0:\) The true mean of the differences is \(\mu_0 = 0\).

We try the following log transformation: \(x' = log(100,100+x)\). Note that we add a constant so that we do not get negative values and then use the base 100 logarithm. The new value of our null hpyothesis is then \(\mu_0' = 1\) because \(x' = log(100,100+x)\) and therefore \(\mu_0' = log(100,100 + \mu_0) = 1\).

library(car)
qqPlot(log(100,100+mouse_cooling$diff))

## [1] 15 18

ggplot(mouse_cooling) + 
  geom_histogram(aes(x = log(100,100+diff)),bins=7)

This looks much better now and we perform a t-test on these data:

results = t.test(log(100,100+mouse_cooling$diff),mu = 1)
results

## 
##  One Sample t-test
## 
## data:  log(100, 100 + mouse_cooling$diff)
## t = -1.4466, df = 18, p-value = 0.1652
## alternative hypothesis: true mean is not equal to 1
## 95 percent confidence interval:
##  0.9254516 1.0137511
## sample estimates:
## mean of x 
## 0.9696013

An alternative would be the sign test:

k = sum(mouse_cooling$diff<0)
results.sign = binom.test(k,n=19,p=0.5)
results.sign

## 
##  Exact binomial test
## 
## data:  k and 19
## number of successes = 7, number of trials = 19, p-value = 0.3593
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.1628859 0.6164221
## sample estimates:
## probability of success 
##              0.3684211

In both cases we find no significant devation from 0 in the differnces. We thus cannot reject the null hypothesis.

Exercise 25

1. \(\chi^2\) Goodness of fit test H0: Fish numbers per unit area follow a Poisson distribution.

2. Two sample t-test or a non-parametric alternative. H0: There is no difference in the mean tumor size change in the two groups.

3. One sample t-test or a non-parametric alternative. H0: Themean tumor size change durign the treatment is 0.

4. \(\chi^2\) Goodness-of-fit test H0: The number of smokers in a city is proportional to the number of people in that city.

5. One-sample t-test H0: Mean change in patient body mass during hospital stay is zero.

Exercise 26

health.expenditure <- read.csv("~/Dropbox/Teaching/StatisticsForBiology/health_expenditure.csv")

ggplot(health.expenditure) +
  geom_histogram(aes(x = Per.capita.health.expenditure.in.2010),bins=10)

Becaue it is highly right skewed.

3. The sample size is

n = dim(health.expenditure)[1]
n

## [1] 18

4. Calculate the natural logarithm of each data point.

data = mutate(health.expenditure,log.expend=  log(Per.capita.health.expenditure.in.2010))
data$log.expend

##  [1] 6.761573 5.768321 5.476464 6.782192 6.156979 4.276666 4.094345 6.408529
##  [9] 3.850148 5.192957 5.509388 4.691348 7.436617 4.997212 5.888878 4.343805
## [17] 7.305860 6.522093

ggplot(data) +
  geom_histogram(aes(x = log.expend),bins=4)

qqPlot(data$log.expend)

## [1] 13  9

5. The mean is

 M =mean(data$log.expend)
M

## [1] 5.636854

6. The standard deviation is

SD = sd(data$log.expend)
SD

## [1] 1.110587

The standard error is

SE = SD/sqrt(n)
SE

## [1] 0.2617679

CI.log = t.test(data$log.expend)$conf.int
CI.log

## [1] 5.084572 6.189136
## attr(,"conf.level")
## [1] 0.95

or alterantively we can approximate it by:

CI.log.app = c(M-2*SE,M+2*SE)
CI.log.app

## [1] 5.113318 6.160390

CI.nonlog = exp(CI.log)
CI.nonlog

## [1] 161.5108 487.4248
## attr(,"conf.level")
## [1] 0.95

M.nonlog = mean(data$Per.capita.health.expenditure.in.2010)
hist(data$Per.capita.health.expenditure.in.2010,col="firebrick",main="",breaks=10)
abline(v = CI.nonlog,lty=2)
abline(v = M.nonlog,lwd=2)

Note how asymmetric the back transformed CI is - this reflects the asymmetry of the data. Be careful however, as this is an approximation! (Note: mean of logarithm transformed data is not the same as the logarithm of the mean of the data).

Exercise 27

no.recycling = c(4,4,4,4,4,4,4,5,8,9,9,9,9,12,12,13,14,14,14,14,15,23)
recycling = c(4,5,8,8,8,9,9,9,12,14,14,15,16,19,23,28,40,43,129,130)

n.rec = length(recycling)
n.no.rec = length(no.recycling)

group = rep(c("Recycling","No.Recycling"),times=c(n.rec,n.no.rec))

dat = data.frame(var = c(recycling,no.recycling),group = group)

ggplot(dat) + 
  geom_histogram(aes(x = var,fill=group))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(dat) + 
  geom_histogram(aes(x = var,fill = group),bins=12) + 
  facet_wrap(~group,nrow=2) 

Both shape and spread seem to be different here, altough it is diffcult to tell based on the small sample sice. A Normal Distribution can however be ruled out, especially fro the recycling group.

ggplot(dat) + 
  geom_qq(aes(sample = var,col=group)) + 
  geom_qq_line(aes(sample = var,col=group))+ 
  facet_wrap(~group,nrow=2) 

qqPlot(recycling)

## [1] 20 19

qqPlot(no.recycling)

## [1] 22 21

2. A t-test should not be used. Based on the small sample size, a permutation test will not work well (see exercise 27). We opt for a Man-Whitney U test.

c and d)

\(H_0:\) It is equally likely that a randomly selected value from one group will be less than or greater than a randomly selected value from a second population. \(H_A:\) It is not equally likely that a randomly selected value from one group will be less than or greater than a randomly selected value from a second population.

wilcox.test(var~group,data=dat)

## Warning in wilcox.test.default(x = c(4, 4, 4, 4, 4, 4, 4, 5, 8, 9, 9, 9, :
## cannot compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  var by group
## W = 126.5, p-value = 0.01825
## alternative hypothesis: true location shift is not equal to 0

# alternative without dataframes
wilcox.test(recycling,no.recycling)

## Warning in wilcox.test.default(recycling, no.recycling): cannot compute exact p-
## value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  recycling and no.recycling
## W = 313.5, p-value = 0.01825
## alternative hypothesis: true location shift is not equal to 0

For comparison, a t-test in this case would yield (not taht the t-test uses a different null hypothesis ):

t.test(var~group,data=dat)

## 
##  Welch Two Sample t-test
## 
## data:  var by group
## t = -2.1447, df = 19.681, p-value = 0.04465
## alternative hypothesis: true difference in means between group No.Recycling and group Recycling is not equal to 0
## 95 percent confidence interval:
##  -34.9245559  -0.4663532
## sample estimates:
## mean in group No.Recycling    mean in group Recycling 
##                   9.454545                  27.150000

Both test actually give similar results, which is a good sign. We therefore reject the null hypotesis. A randomly selected bin from the no recycling group is more likely to contain more paper as compared to a randomly selected bin from the other group.

Exercise 28

We first calculate the actual difference of the means:

no.recycling = c(4,4,4,4,4,4,4,5,8,9,9,9,9,12,12,13,14,14,14,14,15,23)
recycling = c(4,5,8,8,8,9,9,9,12,14,14,15,16,19,23,28,40,43,129,130)

n.rec = length(recycling)
n.no.rec = length(no.recycling)

group = rep(c("Recycling","No.Recycling"),times=c(n.rec,n.no.rec))

dat = data.frame(var = c(recycling,no.recycling),group = group)
MeanR = mean(no.recycling)
MeanNR = mean(recycling)
diffMeans =  MeanNR - MeanR

nPerm = 1000
permResult <- vector() # initializes
for(i in 1:nPerm){
    # step 1: permute the group labels
    permSample <- sample(dat$group, replace = FALSE)
    
    # step 2: calculate difference betweeen means
    permMeanR <-mean(dat$var[permSample=="Recycling"])
    permMeanNR <-mean(dat$var[permSample=="No.Recycling"])
    permResult[i] <- permMeanNR - permMeanR
}

Plot null distribution based on the permuted differences.

hist(permResult, right = FALSE, breaks = 100,col="darkred")

Note how the distirbuion looks odd. This is due to the small sample size: \(n_1 = 20, n_2 =22\). More precily, this is problematic because there are few extreme values in the data set (129 and 130 in the recycling dataset). These exterem points will ahve a strong impact on the empirical null distirbution: if both are in the same resampled group, they will shift the whole distirbution strongly to the left or right. This causes the bimodal distribution. This is a sign that the sample size might be too small for a permutation test. However, usually this leads to small power and large p-values. So we cannevertheless use the null distribution to calculate an approximate P-value, we just need to be aware that our results will be very conservative. The p value is twice the proportion of the permuted means that fall below the observed difference in means, diffMeans. The following code calculates the number of permuted means falling below diffMeans.

sum(as.numeric(permResult >= diffMeans))

## [1] 1

These commands obtain the fraction of permuted means falling below diffMeans.

sum(as.numeric(permResult >= diffMeans)) / nPerm

## [1] 0.001

Finally, multiply by 2 to get the P-value for a two-sided test.

2 * ( sum(as.numeric(permResult > diffMeans)) / nPerm )

## [1] 0.002

We can illustrate the permutaed differences and add the observed difference on to show the outcome of our test: the birght red colored bars indicate the proportion of permutation results that yielded a more exterme results as the one observed.

hist(permResult[permResult<=diffMeans],breaks=seq(-20,20,by=1),col="darkred",main="")
hist(permResult[permResult>diffMeans],breaks=seq(-20,20,by=1),add=T,col="red")
abline(v = diffMeans,col="red",lwd=2)

An alternative approach is

library(perm)
permTS(dat$var[dat$group=="Recycling"],dat$var[dat$group=="No.Recycling"],exact = T)

## 
##  Exact Permutation Test Estimated by Monte Carlo
## 
## data:  GROUP 1 and GROUP 2
## p-value = 0.004
## alternative hypothesis: true mean GROUP 1 - mean GROUP 2 is not equal to 0
## sample estimates:
## mean GROUP 1 - mean GROUP 2 
##                    17.69545 
## 
## p-value estimated from 999 Monte Carlo replications
## 99 percent confidence interval on p-value:
##  1.003509e-05 1.482735e-02

We again find that we can reject the Null Hypothesis.

Exercise 29

dom = "http://www.zoology.ubc.ca/" 
path = "~schluter/WhitlockSchluter/wp-content/data/" 
dat = "chapter15/chap15e1KneesWhoSayNight.csv" 
location = paste(dom,path,dat,sep="") 
circadian <- read.csv(url(location),stringsAsFactors = FALSE) 

1. Visualize the data:

ggplot(data = circadian) + 
  geom_violin(aes(x = treatment,
                  y =shift),
              fill=rgb(0.3,0.5,0.6,0.7),
              bw = 0.5) + 
  geom_point(aes(x = treatment,y = shift),
             col="darkslategray",size=3) +
  theme_minimal()

2. Test for differences in means:

mod = lm(shift ~ treatment,data = circadian)

anova(mod)

## Analysis of Variance Table
## 
## Response: shift
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## treatment  2 7.2245  3.6122  7.2894 0.004472 **
## Residuals 19 9.4153  0.4955                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We find that not all groups have the same mean.

pairwise.t.test(circadian$shift,circadian$treatment)

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  circadian$shift and circadian$treatment 
## 
##      control eyes  
## eyes 0.0088  -     
## knee 0.9418  0.0088
## 
## P value adjustment method: holm

We find that the the “eyes” is different from the other two groups, but “knee” is the same as “control”. This can also be seen very well on the plot of the raw data.

4. The new study failed at replicating the older results, which suggested that light treatment behind the knee can alter your circadian rhythm of melatonin production. Light treatment behind the knee indeed has a different effect than light treatment to the eyes. However, inclusion of a control group shows that only light treatment at the eyes showed an effect relative to the control group. This illustrates how important a proper control group is in experimental design. Without a proper control group, it is nearly impossible to determine whether a certain treatment indeed has an effect.

Exercise 30

Load the data set:

dom = "http://www.zoology.ubc.ca/" 
path = "~schluter/WhitlockSchluter/wp-content/data/" 
dat = "chapter17/chap17e1LionNoses.csv" 
location = paste(dom,path,dat,sep="") 
lion <- read.csv(url(location)) 
head(lion) 

##   proportionBlack ageInYears
## 1            0.21        1.1
## 2            0.14        1.5
## 3            0.11        1.9
## 4            0.13        2.2
## 5            0.12        2.6
## 6            0.13        3.2

ggplot(data = lion) + 
  geom_point(aes(x = ageInYears,y = proportionBlack)) +
  theme_classic()

2. The correlation coefficient can be calcualted using the function cor:

cor(lion$ageInYears,lion$proportionBlack)

## [1] 0.7898272

model = lm( ageInYears~ proportionBlack,data=lion)
summary(model)

## 
## Call:
## lm(formula = ageInYears ~ proportionBlack, data = lion)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5449 -1.1117 -0.5285  0.9635  4.3421 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.8790     0.5688   1.545    0.133    
## proportionBlack  10.6471     1.5095   7.053 7.68e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.669 on 30 degrees of freedom
## Multiple R-squared:  0.6238, Adjusted R-squared:  0.6113 
## F-statistic: 49.75 on 1 and 30 DF,  p-value: 7.677e-08

The coefficients of our model are:

print(model$coefficients)

##     (Intercept) proportionBlack 
##       0.8790062      10.6471194

The linear model for the regression line is therefore

\[ageInYears = 0.879 + 10.647 * proportionBlack\]

library(car)
qqPlot(model$residuals,dist="norm")

## [1] 30 27

A QQ-plot of the residuals shows that they are normally distributed.

plot(model$fitted,model$residuals)

The Tukey-Anscombe plot shows no (strong) deviation from uniform variance of residuals.

We get confindence intervals with the function confint:

confint(model)

##                      2.5 %    97.5 %
## (Intercept)     -0.2826733  2.040686
## proportionBlack  7.5643082 13.729931

This can be directly seen from the CIs from (e) or the summary of the model:

summary(model)

## 
## Call:
## lm(formula = ageInYears ~ proportionBlack, data = lion)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5449 -1.1117 -0.5285  0.9635  4.3421 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.8790     0.5688   1.545    0.133    
## proportionBlack  10.6471     1.5095   7.053 7.68e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.669 on 30 degrees of freedom
## Multiple R-squared:  0.6238, Adjusted R-squared:  0.6113 
## F-statistic: 49.75 on 1 and 30 DF,  p-value: 7.677e-08

The p-value for the null hypothesis that the slope is 0 is very small, \(P = 7.7 * 10^-8\). We can therfore reject the null hypothesis that \(beta_1 = 0\).

This is easiest using geom_smooth() in ggplot:

ggplot(data = lion,aes(x = ageInYears,y = proportionBlack)) + 
  geom_point() +
  geom_smooth(method = "lm") +
  theme_classic()

## `geom_smooth()` using formula 'y ~ x'

See what happens when we do not specify the method in lm():

ggplot(data = lion,aes(x = ageInYears,y = proportionBlack)) + 
  geom_point() +
  geom_smooth() +
  theme_classic()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'